function [mom_ave] = simulate(ave_theta, h_theta, h_eps, delta, cost,gamma, wc, SIM,step,learn)
% solve for each possible theta draw
global num_firm T_sim S T n_theta  beta n_tran tau mktfe search T_cal country_num phi 

%% solve the optimal price path
% discretize possible firm draws
sigma = 0.17* cost*exp(ave_theta);
theta_firm = linspace(ave_theta-2.5*sqrt(1/h_theta)-2.5*sqrt(1/h_eps),ave_theta+2.5*sqrt(1/h_theta)+2.5*sqrt(1/h_eps),n_theta);
theta_firm_f = [theta_firm(2:n_theta) theta_firm(n_theta)];
theta_firm_l = [theta_firm(1) theta_firm(1:n_theta-1)];
pr_theta_firm = normcdf((0.5*(theta_firm+theta_firm_f)-ave_theta)*sqrt(h_theta))-...
    normcdf((0.5*(theta_firm+theta_firm_l)-ave_theta)*sqrt(h_theta));
pr_theta_firm = pr_theta_firm/sum(pr_theta_firm); %1*n_theta, determine which true theta the firm draws

price_policy = nan(n_tran*country_num,T+1,n_theta);
theta_grid_mat = nan(n_tran,T+1,n_theta);
pr1_mat = nan( n_tran*country_num,n_theta);
pr2_mat = nan(n_tran*T, n_tran*country_num, n_theta);
for i = 1:length(theta_firm)
    [p theta_grid pr1 pr2] = solvemain(ave_theta, theta_firm(i), h_theta, h_eps, delta,cost,wc);
    if sum(sum(isnan(p)))>0
        sprintf('%2.4f, %2.4f,%2.4f,%2.4f,%2.4f,%2.4f',ave_theta, theta_firm(i), h_theta, h_eps, delta,cost)
        sprintf('%2.2f theta_n gives nan price policy\n',i)
    end
    price_policy(:,:,i) = p;% n_tran * country_num by T+1; this is price with delivery fee
    theta_grid_mat(:,:,i) = theta_grid;% n_tran by T+1 by nfirm;for each firm
    pr1_mat(:,i) = pr1'; % column is different firm, row is first period theta state after country impression, n*country_num by nfirm
    pr2_mat(:,:,i) = pr2; % row is from zero period to T period for all possible theta state; column is next period's state, n*T by n*country_num
end

%sprintf('DP problem solved---------------------------------------\n')

%replace_tmp = nan(size(SIM.replace));
exit_tmp = nan(size(SIM.exit));
for s = 1:S
    exit_tmp(:,:,s) = transform(repmat([delta;1-delta],num_firm),SIM.exit(:,:,s)'); % 1 = exit; 2= stay
end
%SIM.replace = replace_tmp;
SIM.exit = exit_tmp;

%% simulation

% - dead firm replacement matrix, num_firm by T_sim
% - consumer country matrix
% - firm exit matrix
% - consumer buying matrix
% - firm ads matrix
% - firm feedback matrix
mom_all = [];
for k = 1:S
    SIMINPUT = struct('sim_matrix_num',6);
    SIMINPUT.replace = SIM.replace(:,:,k);
    SIMINPUT.country = SIM.country(:,:,k);
    SIMINPUT.exit = SIM.exit(:,:,k);
    SIMINPUT.purchase = SIM.purchase(:,:,k);
    SIMINPUT.ads = SIM.ads(:,:,k);
    SIMINPUT.feedback = SIM.feedback(:,:,k);
    if learn==1
        [mpurchase mprice mom sample_country sample_net_sale ann_sale ann_price ann_purchase ann_rating] = simulate_one(SIMINPUT,T_cal, delta, h_theta, h_eps, ave_theta, ...
        theta_firm, theta_grid_mat, pr_theta_firm,pr1_mat,pr2_mat, price_policy, cost, gamma,step);
    else
        [mpurchase mprice mom sample_country sample_rating ann_sale ann_price ann_purchase ann_rating] = simulate_one_nolearn(SIMINPUT,T_cal, delta, h_theta, h_eps, ave_theta, ...
        theta_firm, theta_grid_mat, pr_theta_firm,pr1_mat,pr2_mat, price_policy, cost, gamma,step);
    end
    mom_all = [mom_all mom];
     filename1 = sprintf('month_price_#%d.txt',k);
     filename2 = sprintf('month_sale_#%d.txt',k);
      filename3 = sprintf('month_rating_#%d.txt',k);
     dlmwrite(filename1,ann_price);
     dlmwrite(filename2,ann_sale);
     dlmwrite(filename3,ann_rating);
 end
mom_ave = nanmean(mom_all,2);
%sprintf('%2.2f simulations done!-----------------------------\n',S)
end
